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We analyze correlations in step-edge fluctuations using the Bortz-Kalos-Lebowitz kinetic Monte 
Carlo algorithm, with a 2-parameter expression for energy barriers, and compare with our VT- 
STM line-scan experiments on spiral steps on Pb(lll). The scaling of the correlation times gives 
a dynamic exponent confirming the expected step-edge-diffusion rate-limiting kinetics both in the 
MC and in the experiments. We both calculate and measure the temperature dependence of (mass) 
transport properties via the characteristic hopping times and deduce therefrom the notoriously- 
elusive effective energy barrier for the edge fluctuations. With a careful analysis we point out the 
necessity of a more complex model to mimic the kinetics of a Pb(lll) surface for certain parameter 
ranges. 
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I. INTRODUCTION 

Below the roughening temperature a rich variety of 
structures are observed on crystal surfaces: islands, 
mounds, shallow and deep pits, self-assembled patterns, 
step bunches, and spiral steps, as well as such extended 
structures as quantum wires and arrays of quantum dots. 
The evolution of such nanostructures proceeds via mass 
transport on flat terraces and along the steps separat- 
ing themi^ The motion of these structures can now be 
directly observed and quantitatively determined by mon- 
itoring the constituent steps using advanced experimen- 
tal techniques like variable temperature scanning tun- 
neling microscopy (VT-STM) with near-atomic spatial 
resolution and with a time resolution fast enough to ob- 
serve fluctuations on the nanometer scale (but too slow 
to observe individual atomic movements, as we shall sec 
below). The kinetic parameters relevant to the mass 
transport mechanisms are typically determined from the 
spatio-temporal evolution and/or fluctuation of steps. 
Thus, scrutiny of the relatively simple behavior of steps 
offers a gateway to understanding the evolution of com- 
plex structures and an avenue to test theoretical argu- 
ments against the (harsh) reality of experiments. 

In this paper we use STM to study thermal fluctuations 
of spiral steps^ on the surface of Pb(lll) crystallites A 
The principal mass transport mechanism on these sur- 
faces is known to be periphery diffusion^ near and mod- 
erately above room temperature. Less is known, however, 
about the temperature dependence of the mass transport, 
with which an effective energy barrier can be associated 
in an Arrhcnius-type picture due to the activated nature 
of these processes. This effective energy barrier, which 
can be measured in experiments, is intimately related to 
more fundamental parameters characterizing the surface 
and its steps: the kink energy and the energy barrier in 
diffusion. These parameters describe both the equilib- 



rium and out-of-equilibrium behavior of the system, e.g. 
the equilibrium crystal shape, and the linear response of 
the system to small perturbations, the diffusion coeffi- 
cient. 

We set up a theoretical model and subject it to exten- 
sive kinetic Monte Carlo (kMC) simulations. In contrast 
to a step on a vicinal surface, the azimuthal orientation^ 
of a spiral step is not well-defined but instead spans all 
angles. In the experiments an average is taken over the 
undefined directions in an imprecise way. In our model 
we take a simple average of all azimuthal directions, as 
one would for an island^ These orientation-averaged fluc- 
tuations provide the basis for calculating the transport 
properties and the associated effective energy barriers in 
this paper. 

The paper proceeds as follows. We first recall the im- 
portant results of Langevin analysis applied to step fluc- 
tuations dominated by atomic motion along the steps. In 
the third section we describe the experimental technique 
for measuring spirals on Pb crystallites and our results. 
In Section IV we describe our extensive kMC simulations 
on a hexagonal lattice using energy barriers determined 
by a subset of occupied nearest-neighbor sites of each 
hopping atom with energies computed with a semiempiri- 
cal scheme optimized for surfaces. (An appendix expands 
on the issue of detailed balance in such models.) Section 
V presents our kMC results and comparisons with exper- 
imental values, followed by a final section of discussion 
and conclusions. 



II. BACKGROUND THEORY 

Step-edge fluctuations are typically described with 
stochastic partial differential equations drawn from 
Langevin formalism. The characteristic parameters are 
equilibrium values determined phenomcnologically or 
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from simple models such as a lattice-gas (Ising-likc) 
model of the appropriate symmetry. We detail here 
the necessary equations that apply to the fluctuations 
of steps on the Pb(lll) surface. As we will see, the fluc- 
tuations are due to step edge diffusion, which is known 
to be described by the Cahn-Hilliard equation 
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where x = x(y, t) is the perpendicular displacement of 
the step at point y along the step at time t, Th is 
the step mobility, which encompasses the friction felt 
by the step when randomly kicked, ft is the step stiff- 
ness, the inertial parameter which counterweights devi- 
ations of the step from its equilibrium shape. The con- 
served thermal noise r]h{y, t) satisfies (r)h{y, t)r]h{y', t')) = 
—2Tf l V 2 S(y — y')5(t — t'), where () is the canonical ensem- 
ble average. Here we neglect the fact that the observed 
objects (the spiral and island) have an intrinsic curva- 
ture, but take them as locally straight and apply the 
same formula as for the straight step. This amounts to 
leaving out a geometric factor in Eq. that is close to 
unity. 

The time correlation function G(t) = ([x(Q) — x(t)] 2 ) 
is readily obtained^ from the Cahn-Hilliard equation for 
times shorter than the correlation time of the system 
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where the Euler gamma function T(3/4) w 1.2254. We 
denote by c(T) the prefactor of t 1 / 4 ; it is the number we 
typically extract as the endproduct of our experimental 
measurements. 

The step mobility is directly related to the microscopic 
time constant Th that indicates the average time between 
moves at any site (similar to the collision time in trans- 
port theory) . An equation similar to the Nernst-Einstein 
relation for the diffusion coefficient describes this rela- 
tionship = £l 5 / 2 /i7j, where il is the area of the 2d 
unit cell a Pb atom occupies in the (111) surface. Thus, 
in experimental analyses based on Eq. @, the time con- 
stant is 
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Though the above equations are all continuum equa- 
tions, the step stiffness (3 is typically drawn from a simple 
lattice-gas model. The transition from continuum to lat- 
tice or vice versa implicitly contains the assumption that 
the steps in question are long enough to be able to form 
capillary waves with negligible finite-size effects. We cal- 
culated the stiffness from the Akutsu formula-^ for (111) 
surfaces or its low temperature approximation^ 



Here a is the lattice constant of the hexagonal lattice, 
and £fc is the kink energy. In the present study we con- 
sider only nearest-neighbor (NN) interactions between 
atoms (in particular when computing hopping barriers) 
with [attractive] interaction energy —Enn, i.e. with an 
effective bond energy Enn > 0. In this basic lattice-gas 
model, there is the simple relation e k = Enn/2 between 
these characteristic energies^. 

The STM measurements are always so-called line-scan 
measurements: the fluctuations are measured only along 
a single line perpendicular to the step edge. Since this 
line is fixed, we cannot gather information about a longer 
stretch of the edge as one would typically in LEEM mea- 
surements and in Monte Carlo simulations. While this 
low information density is a severe limitation of these 
measurements, it is amply compensated by the very high, 
basically atomic resolution, in contrast to LEEM mea- 
surements. 

In the MC simulations, since we have data for the 
whole spatial extension of the step (or island in our 
case) we can analyze the (spatial) Fourier spectrum of 
the fluctuations. In this case x(y, t) measures the de- 
viations from the equilibrium 2d crystal shape due to 
the fluctuations and x(y,t) — ^ q x q {t) exp[igy], where 
q = u2t:/L, L is the circumference of the island; n = 
-N/2 + 1, -N/2 + 2, N/2, where N is the number of 
degrees of freedom of the island edge, i.e. the number of 
perimeter atoms. The time correlation function of the 
Fourier modest G q (t) = (\x q (t) — x q (Q)\ 2 ) exponentially 
approaches the equilibrium fluctuation (roughness) of the 
step 
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The correlation time r(q) scales with the wave number 
as 

k B T T h 
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with the dynamical exponent z = 4 in the periphery dif- 
fusion case, where Th is the same hopping time as above. 

Since the kinetic processes are activated processes on 
the surface, the hopping rate 
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= (2kT/3a)exp[e k /kT]. 
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can be characterized by an effective hopping barrier Eh 
and an attempt frequency v. The hopping barrier com- 
bines the kink energy e k responsible for the kink density 
along the step and the (edge-)diffusion barrier Ed- For 
any particular crystallographic direction, this combina- 
tion comes with definite prcfactors. However, because we 
average over all directions, this decomposition becomes 
irrelevant; therefore, we directly measure the hopping 
barrier both in the experiment and the MC simulations. 
We expect the attempt frequency to be of order the De- 
bye frequency since it is mostly the lattice vibrations that 
give rise to adatom hops on the surface. 
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FIG. 1: (Color online) lOOOnm x lOOOnm STM image of a 
spiral dislocation on top of a Pb(lll) crystallite. The STM 
line-scan measurements were performed perpendicular to the 
step edge, typically 5-50nm from the dislocation core. 



III. EXPERIMENTAL METHODS AND 
RESULTS 

Supported 3D Pb crystallites were prepared on 
Ru(0001) by dewetting a continuous film (10-30 nm nom- 
inal thickness) under UHV conditions. The film was pre- 
pared by depositing Pb from a heated alumina tube onto 
a room-temperature Ru(0001) substrate. By heating the 
sample to the melting point of Pb (600 K) for several 
minutes, the continuous film formed liquid beads that 
were then quenched to form the solid crystallites. Crys- 
tallites formed in this way have (111) facets parallel to 
the Ru substrate plane that can be as large as 1 /Ltm in 
diameter. 

STM measurements were made with a commercial in- 
strument (Omicron VT-STM) at room temperature and 
above. The sample was heated by a PBN heater in the 
sample holder (Omicron), and the temperature was cal- 
ibrated using a manufacturer-supplied curve that was 
checked with infrared pyrometry for T > 520 K. Step 
fluctuations were monitored on (well isolated) steps us- 
ing the line-scan method of imaging: the STM tip is fixed 
over one point along the step edge and scanned perpen- 
dicularly for 20 - 120 seconds. The step position from 
such an image could be extracted and used to compute 
G(t), as described in Section II. 

We make our measurements in the 300K-400K temper- 
ature range. The experimental situation strongly changes 
in this range. At high temperatures (> 350K), kinet- 
ics are fast and small structures accidentally created on 
the facet such as islands, voids or mounds, decay away 
fairly fast. Isolated steps are only obtained when there 



is a screw dislocation on the facet as shown in Fig. ^ 
The spiral step terminating at the screw dislocation con- 
tinues to unwind until it reaches equilibrium with the 
crystallite. For temperatures < 350K, this unwinding 
process becomes slow, so that it is more difficult to de- 
termine how close the step is to equilibrium. As for room 
temperature, small structures do not decay easily. Thus 
for room temperature experiments, in addition to spiral 
steps, we consider steps crossing in-between small struc- 
tures on the facet, treating them as isolated steps which 
are assumed to be in local equilibrium. 

Wc have made two sets of measurements, one with less 
temperature resolution, in which we explore the coarse 
structure of the temperature dependence, and another 
with higher resolution, in which wc pay more attention 
to the details. We applied different analyses to the two 
sets; there are also differences in the experimental envi- 
ronment. The effects of these differences will be discussed 
later. 

Experimentally determined temporal correlation func- 
tions (from the higher-temperature resolution measure- 
ments) for different temperatures are shown as a log- 
log plot in Fig. [3 The points are the results of aver- 
aging several (5-15) experimental G(t) curves obtained 
from displaccment-vs.-time images taken over typically a 
few minutes at the same point of the step edge for each 
temperature. For screw dislocations most measurements 
were done relatively close to the dislocation core (within 
5-50nm, which is also far enough from the core, so that 
the screw step has reached its full value), although we do 
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FIG. 2: (Color online) Log-log plot of equilibrium correlation 
functions (measured in A 2 ) vs. time (sec). Measured in the 
configuration depicted in Fig. the data sets are for tem- 
peratures 300K, 320K, 330K, 340K, 350K, 370K, 390K, from 
bottom to top, respectively. The slope of the curves give the 
growth exponent and determine the rate limiting process re- 
sponsible for the fluctuation of the steps. The thick dashed 
line corresponds to an exponent 1/4. Evidently the data for 
the temperatures depicted in this plot are essentially consis- 
tent with this value (see text) , suggesting that the rate limit- 
ing process is step-edge diffusion. 



4 



not see a significant difference in the trial measurements 
farther from the core. 

The G(t) data are linear on log-log plots, indicating 
power-law growth of G(t); their slopes are consistent 
with the t 1 ' 4 behavior (indicated by the slope of the 
dashed line in Fig. 0) expected for step-edge diffusion 
(cf. Eq. (J2J). Specifically, their exponent values are close 
to 1/4, ranging from 0.20 to 0.35 and clustering around 
0.27. Consistent with Eq. ©, the magnitude of G(t) (i.e. 
the prefactor c(T)) increases dramatically as temperature 
increases. 

From a fit of the experimental curves like the ones in 
Fig. El using Eq. J5J, c(T) is obtained. In Fig.Olwe show 
the experimental prefactors as a function of temperature. 
Each data point is obtained from an average of all G(i) 
curves for each measured step position. The set of mea- 
surements with little temperature resolution is indicated 
by triangles, while the ones giving more detail about the 
temperature dependence are denoted by squares. The 
latter set consistently gives somewhat higher prefactor 
values probably due to the different image processing and 
analysis used. 

We use the prefactors to extract the step-edge diffu- 
sion time constant Th with the aid of Eq. [3J It is found 
to decrease monotonically from 50/is to 0.073ns as tem- 
perature increases from 300K to 400K. Due to the fourth 
power of c(T) in Eq. (jSJ), the error in the determination 
of Th can be very significant. In later sections this exper- 
imental uncertainty will be discussed more fully. These 
findings for Th are discussed in Sec. |V] and depicted in 
Fig. □ along with the MC results. 
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FIG. 3: (Color online) Prefactors c(T) of the experimentally 
measured correlation functions vs. temperature, calculated 
based on Eq. [5] assuming step-edge diffusion as the rate lim- 
iting process for the fluctuations of the steps. For each tem- 
perature, the different data points correspond to the averaged 
prefactor at each measured step position. The triangles and 
squares indicate two different sets of measurements (see text). 



IV. KINETIC MONTE CARLO 

In order to investigate correlations in a "clean" envi- 
ronment, we perform a kinetic Monte Carlo simulation 
of a hexagonal lattice gas. Simulational data can pro- 
vide corroboration of our fundamental picture of what 
transpires during the experiments and can also help in 
their analysis. It has the advantage that far more data 
can be collected, allowing for an even more precise anal- 
ysis, and it also has spatial resolution {along step edges) 
which would be difficult to get in today's STM measure- 
ments, as will become more evident shortly. While these 
measurements can acquire fluctuation data in a line scan, 
which provides us with fluctuation in time at a particu- 
lar point along the step edge, the kMC provides us with 
instantaneous snapshots of a spatially extended region of 
the step edge (again along the step) , in our case snapshots 
of the whole island edge. Data similar to what we can get 
in the kMC simulation can be obtained by LEEMji^ but 
in that case the spatial resolution both along and per- 
pendicular to the step is less than the resolution of the 
simulation (or the resolution of the STM perpendicular 
to the step edge). Thus, both experimental methods have 
their advantages and disadvantages over the other while 
the advantages of the MC mentioned above is diminished 
by the fact that it inevitably contains simplifications in 
the geometry of the system and the energy landscape of 
the surface determining the kinetic barriers for the hops 
of the adatoms and also for the configuration of the hops. 
(Only NN hops of single adatoms are considered in the 
simulations.) 

In considering the (111) surface of Pb to be a 2D, 
hexagonal net with 6-fold symmetry, we disregard the 
ABC stacking structure leading to only 3-fold surface 
symmetry. The adatoms on the surface are put into a cir- 
cular "container" with an impermeable wall. The adatom 
gas in the container is a highly supersaturated (but 
low-density) lattice gas at the temperatures of present 
concern; hence, most of the atoms in the 2D gas con- 
dense into a big island close to the middle of the con- 
tainer, surrounded by detached adatoms moving around 
the terrace between the island and the container wall. 
This very rarified remaining gas of mobile adatoms and 
edgcatoms (adatoms sliding along the island edge) gives 
rise to the fluctuations of the island edge. The island is of 
monatomic height; we completely neglect the possibility 
of an adatom hopping onto the top of the island. 

We use the Bortz-Kalos-Lebowitaii (BKL) continuous- 
time MC algorithm in our simulations since it is well 
suited to the low temperatures of interest here. The crit- 
ical temperature in the simulations is T c = 1400K, while 
the typical temperatures in question are below 400K. The 
rejection-free feature of the BKL algorithm greatly im- 
proves the efficiency of the simulations, and efficiency is 
further enhanced by the n-fold-way method, since n in 
our case is quite low (n = 5) due to the small number of 
different energy barriers. 

We evaluate energy barriers using scmicmpirical sur- 
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face embedded atom method^ (SEAM) calculations us- 
ing Ercolessi's glue potentials^ for the Pb(lll) surface 
(fourth column of Tabic |IJ . Our BKL simulation of the 
lattice gas expresses these barriers in terms of two char- 
acteristic energies. One of these parameters is the bar- 
rier for an adatom to hop on the terrace when it has no 
lateral neighbors, i.e., surface diffusion; a direct compu- 
tation of this energy barrier with SEAM gives Etd=7Q 
meV (see Table QJ. The other parameter accounts for the 
possibility that the adatom has lateral neighbors with 
which "bonds" [direct and indirect] are broken when the 
adatom hops. We consider only NN bonds as depicted in 
Fig.H 

In our first parametrization, we adopt a scheme we call 
"break-3" in which the energy barrier of a hop depends 
on how many sites are occupied among the sites "behind" 
a hop, labelled 0, 1, or 1'; they each contribute i?NN=T30 
meV, which is the other characteristic barrier of the lat- 
tice gas. -Enn is obtained via trial and error comparison 
with EAM barriers; see Table |U and the next paragraph. 
We neglect any effect of the other nearest-neighbor sites 
to the two sites involved in the hop (2, 2', 3, 3', and 4), 
which we expect to be smaller. The actual energy bar- 
rier Ef, for such a hop is the sum of the terrace diffusion 
barrier and the contributions of the bond breakings 

E b =E TD +mE NN , (8) 

where m is the number of occupied "behind" sites with 
which a bond is broken, and it ranges from to 3. 

On closer inspection the energy barriers computed as 
above with SEAM (or any form of EAM) seem not to sat- 
isfy detailed balance for more complex processes such as 
the one depicted in Fig. [5] and analyzed in the Appendix. 
The reason is that these processes involve next-nearest- 
neighbor (NNN) and possible further-neighbor (or multi- 
neighbor—) interactions; thus, EAM only satisfies de- 
tailed balance if these are also taken into account, so 
that the hopping process is considered in its full com- 
plexity. However, to allow for long kMC runs and trans- 




FIG. 4: (Color online) Sample configuration on triangular 
lattice, indicating our nomenclature for the nearest-neighbor 
sites of the adatom before and after the move. In this case 
the move is labelled 01 since site and 1 are occupied among 
the eight NN sites surrounding the two sites involved in the 
hop. See Table □ for the energy barriers of the various hop 
configurations. 
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TABLE I: Tabulation of energy barriers (all in meV) for 
selected hopping processes. See Fig. [I] and text for explana- 
tion of the configurations. The 3 bond-breaking (bb3) scheme 
barriers, 5 bond-breaking scheme (bb5), the EAM energy 
barriers, the EAM energy difference of two states, and the 
Kawasaki-type barriers are listed, respectively. We use the 
bb3 scheme in the kMC simulations in this paper. 



parent results, we only want to capture the main hopping 
mechanisms relevant to the material; the full details of 
the energy landscape are an uninteresting complication 
and distraction. We use the EAM barriers as a guide 
to deduce Etd and Enn and in this way simplify them 
into classes such as in the break-3 scheme just described 
(column two). This scheme along with the other simpli- 
fied schemes described below satisfies detailed balance as 
shown in the Appendix. 

Another such classification is the break-five scheme in 
column three, which involves sites 0, 1, 1', 2, and 2'. 
Each of these sites, if occupied, contributes Enn to the 
hopping barrier, so that in Eq. © m ranges from to 
5. We expect that the break-3 and break-5 parameteri- 
zations bracket the actual dependence of the barrier on 
lateral occupation. 

The break-3 classification is evidently quite a good ap- 
proximation to the EAM in the 0, 12, 012, 011', 23, 11' 
cases. The break-5 is good in the 0, 011', 11' and espe- 
cially the 011'2 cases but in several other cases signifi- 
cantly worse than the break-3 scheme, e.g. 123, 1234, 23, 
and 22' cases. 

For the sake of comparison, we also list the Kawasaki- 
type energy barriers^ (column six) that are usually cal- 
culated based on how many bonds have been broken and 
how many new ones have been formed in a move, plus 
the usual terrace diffusion barrier. This scheme is often 
used in surface simulations. Its philosophy is markedly 
different from the previous two bond-breaking schemes 
because it takes into account the newly formed bonds, 
i.e. the energy of the final state. The Kawasaki scheme is 
not regarded the best scheme for surface simulations ^ 
but in our case Table [I] shows it to be quite close to the 
EAM. 
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V. RESULTS 

To check the precision of the kMC procedure, we com- 
pared the equilibrium island shapes to analytical results 
by Zia— for a triangular latticed In Fig.[S]these shapes 
are depicted for T=250K and show excellent agreement. 
The slight difference between the MC shapes and the ex- 
act shape may be due to the fact that the shapes are 
always measured from the instantaneous center of mass, 
which may give rise to more rounded shapes in the MC. 
Furthermore, since the islands are not large compared to 
the lattice spacing a, there might well be finite size ef- 
fects (which we have not attempted to analyze here) , but 
this effect should also appear in the experimental stud- 
ies of nano-islands having the same size as in the MC 
simulations. 

From the MC data we calculate the line tension and 
the step stiffness in the two principal directions of the 6- 
fold symmetric triangular lattice making 30° angle. They 
compare favorably (see Table |nj) with exact results and 
with measured values (of line tension) for Pb(lll);2I in- 
dicating that the Ising model and the corresponding lat- 
tice gas model with the chosen energy barriers describe 
these structures well. The deviations of the KMC from 
the exact data and the measured ones basically show the 
error bars and the limit of applicability of the model and 
of the simulations. 

In our KMC simulations the step fluctuations are 
driven by step-edge diffusion throughout the studied tem- 
perature range. The scaling of the correlation time r(q) 
with wavenumber gives a dynamic exponent z = 4 char- 
acteristic of these processes (see Fig. [HJ ■ At higher tem- 



40.5 1 ■ . 1 . ■ 1 . . 1 . r 




38 5 1 ' 1 1 1 ' 1 1 1 1 1 1 

71/12 TC/6 7C/4 

FIG. 5: (Color online) The r vs. <f> plot of the equilibrium 
crystal shape at T=250K, for radii r=20a, r=40a, and the 
exact shape in the continuum limit, denoted by dashed (red), 
thin solid (black), and thick solid (blue) lines, respectively. 
The graphs are appropriately rescaled so that the islands have 
the same size to make shape comparison possible for islands 
of different physical sizes. (This transformation is a mathe- 
matical similarity transformation which does not change the 
shape of the islands, only their size.) The inset shows the 
minimum of the curves magnified for better comparison. 
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TABLE II: Tabulation of the line tension and the rela- 
tive stiffness to test the KMC simulation on the exact lat- 
tice gas results by Akutsu and Akutsui--» In comparison, 
experiment finds 29.65meV/A for 323K, 28.65meV/A for 
353K, 27.2meV/A for 393K for the line tension (i.e., the av- 
erage of that of A and B-type steps, since we neglect their 
difference in the simulations). 

peratures we might expect a crossover to other pro- 
cesses like terrace diffusion or attachment-detachment 
with z = 3 or 2 dynamic exponents, respectively. These 
processes could show up due to the increased roughness 
of steps at higher temperatures, which give rise to sites 
in the step edge that are bound to the step by only two 
bonds instead of the more typical three or four, driving 
the system into the attachment-detachment rate limit- 
ing regime (see RehjUHH). But this scenario is con- 
firmed neither by the experiments^ nor the MC simula- 
tions here. 




FIG. 6: (Color online) Dynamical scaling for the KMC data. 
Scaling of the relaxation time r n versus the wavenumber n 
for T =250K (circle), 300K (square), 350K (diamond), 400K 
(triangle). For the lowest temperature r = 40a, in the other 
cases r = 20a. In all four cases the dynamical exponents, 
given by the slopes of the curves, suggest step-edge diffusion 
limited kinetics, see Table ITTT1 The straight dotted-dash and 
the dashed lines are guides (at arbitrary ordinate) to show 
the expected slopes for step-edge diffusion (z=4) and terrace 
diffusion (z=3), respectively. 

The hopping time characterizes the average time 
that elapses between successive "events" at a site, viz. 
the arrival or departure of an adatom. This change of 
occupation can be part of either a step-edge diffusion pro- 



cess (which is the more typical here) or an attachment- 
detachment process associated with terrace transport. 
Using the exact step stiffnesses and the correlation times 
r(q) of the shortest wavelengths that we measured in the 
MC simulations in Fig. we calculate the more funda- 
mental time constant th for each temperature from the 
correlation time of the Fourier modes via Eq. 

As shown in Fig. [7| the Arrhenius plot of the Monte 
Carlo Th gives the effective energy barrier Eh— 380 meV 
and a frequency vq = 7.0 • 10 12 Hz, which is in the range 
of the Debye frequency. We also plot the two sets of 
experimental r/j, in both cases using two different kink 
energies, 50 meV and 65 meV, in j3 in Eq.|H Here 50 meV 
is the average measured kink energy of the A- and B-typc 
steps, while 65 meV is the kink energy we get in our EAM 
calculations. The two experimental data sets with both 
kink energies basically provide us with the same effective 
energy barrier. The two sets are offset by an order of 
magnitude which gives different prefactors, v. This dif- 
ference might be due to the earlier mentioned different 
analysis used in both sets of measurements and shows 
the susceptibility of the prefactors and the robustness of 
the energy barrier. 

Though the two data sets seemingly give similar re- 
sults for the effective energy barrier, the data set with 
higher temperature resolution reveal more subtleties, so 
we analyze this data set here. For low temperatures the 
experiments apparently bracket the MC with the same 
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TABLE III: Tabulation, for the four temperatures used in the 
KMC simulations, of the dynamical exponent z, system size 
(viz., circumference) L , the wavenumber n and wavelength A 
used in calculating the correlation times, the correlation time 
r(A), and the hopping time r n . Due to the activated nature 
of the processes, the hopping times change by a significant 
factor over the given temperature range. 



energy barrier, while for higher temperatures this exper- 
imental data set tends to bend below the MC, suggesting 
a higher effective barrier. If we crudely fit the experimen- 
tal data without using the MC as a guide, we get 780 meV 
and 735 meV, using the two kink energies, respectively. 
These effective barriers for experiments at higher tem- 
perature seem to differ greatly from the nice scaling of 
the various energies for comparable materials such as Pt, 
Ag and Cu (see Table ITvl and Sec. VI), suggesting the 
conclusion that a phenomenon that is not relevant on 
the other surfaces — also not captured by the (nearest- 
neighbor single-atom-hop) MC algorithm comes into 
play on the Pb(lll) surface. We elaborate on this idea 
in the following Section. 



10 | 1 1 1 1 1 1 1 1 1 r 




l/k B T (eV"') 

FIG. 7: (Color online) Arrhenius plot of hopping times. Two 
sets of experimental results (triangles and squares) are plot- 
ted with two different kink energies, efe=50meV (downward 
pointing triangles, and diamonds) and 65meV (upward point- 
ing triangles, and squares) (see Eq. ©). The KMC data are 
depicted in circles with a thick solid-line as a guide to the eye. 
The error bars are indicated on the experimental data along 
with a linear fit (dashed lines) on the higher temperature reso- 
lution data. The dotted lines indicate how we speculate these 
data might continue if the experiments had gone below room 
temperature. The slopes of the fits (dashed lines) to all data 
points except the highest temperature ones show the energy 
barriers to be apparently higher than that of the MC (see text 
for discussion, and also Table ITvt . 



VI. DISCUSSION AND SUMMARY 

The diffusion barrier along a straight step is usually 
deduced from measurements^ to characterize step-edge 
diffusion kinetics. Its measurement is based on consider- 
able amount of assumptions such as low edge atom den- 
sity, which also entails slightly fluctuating steps. Both 
the fluctuation and the geometric kink density is con- 
sidered to be low. With these assumptions the corre- 
lation function G{t) = (2/7r)r(3/4)Pf /4 (c<rf Arf) 1/4 * 1/4 
where Pk ~ 2 cxp[— ek/ksT] is the kink density, c s t = 
cxp{—E st /kBT] the adatom density at the step edge, 
and D st = Do exp[— E^/ksT] is the tracer diffusion 
along a kinked step edge. A common simple assump- 
tion is that the adatom creation energy at the edge is 
E st = 2efc, and the kinked step diffusion energy barrier 
is related to the straight step-edge diffusion energy bar- 
rier via E\ = 6k + Ed- Do is the prefactor of the dif- 
fusion coefficient. Since the kink energy is usually 
known from other types of measurements, the diffusion 
barrier Ed can be determined by analyzing the G(t) cor- 
relation function since 4E t = 3efc + E st + E% = 6ek + Ed, 
where E t is the measured activation barrier of G(t) oc 
cxp[— Et/kBT]^/ 4 . In Table ITvl these diffusion barriers 
are tabulated for Pt, Cu, and Ag, for which these ener- 
gies are known, along with the kink energy and the bulk 
cohesion energy for the sake of comparison in trends. 

In our case, we extract part of the above G(t), namely 
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TABLE IV: Tabulation of the cohesion energy E co h 
[eV/atom], and for the (111) surface the kink energy [meV], 
and the diffusion barrier along straight steps Ed [meV]. For 
these four materials there is a clear trend both in the exper- 
imental (Exp) and theoretical (Th) values of these energies, 
where Pb proves to be the "softest". For Ed for Pt, Cu, Ag 
the theoretical values are from embedded atom (EAM) and ab 
initio (density functional) calculations. For Ed for Pb we give 
our KMC result, which is corroborated by our experiments at 
lower temperatures; this value excellently follows the trend of 
the other columns of the table. 



Th = l/c s tD st , and in the above approximation our hop- 
ping energy is Eh = 3efe + Ed- From the MC simula- 
tion, for the diffusion energy we deduce ^=185meV , 
which is close to the value Ed=200meV predicted from 
our bond-breaking scheme with our semiempirical input 
energies. The under-10% difference shows how effective 
the assumptions are in the above argument. In spite of 
the fact that in the MC simulation we consider a small 
island with very large curvature and obviously many geo- 
metric kinks, the deduced diffusion barrier along the step 
edge is very close to the prediction based on a straight 
unkinked edge. 

For the lower end of the experimentally measurable 
range (viz. room temperature), our simulational results 
agree with these experiments. It is unfortunate that go- 
ing lower temperatures in the experiments was not possi- 
ble because of the very slow equilibration of the samples. 

In similar analyses of the highcr-tcmpcraturc experi- 
ments, this barrier increases by 400meV. Such a barrier is 
twice as high as for Cu, closer to Pt. Since this is strongly 
counter to the trend in the cohesion energy and the kink 
energy, we should consider other phenomena that are re- 
lated to the special spiral geometry in our measurements 
or ones not present for the other materials. These higher 
energy barriers for higher T may be an indication of dif- 
ferent regions of transport mechanisms on Pb(lll). The 



included in our MC as it requires next nearest neighbor 
hops. 

The effective energy barrier may well depend on the ac- 
tual geometry of the sample. Different curvatures might 
lead to different energy barriers due to change in the 
density of geometric kink sites. It may also depend on 
whether the step edge is compact or more fractal-like 
since the energy barrier contains a weighted combina- 
tion of the few typical bond-breaking barriers where the 
weights depend on the average local geometry around 
hopping adatoms on the step edge. This argument also 
implies that these effective energy barriers may even 
depend on temperature since the local geometry may 
change since the roughness increases with temperature. 
However, these geometric arguments do not seem to be 
corroborated by the MC simulations, as we have just 
seen. 

The larger scatter in the measurement data both in 
this work and also in other's data sets2ii22iiS indicate the 
difficulty of measuring these barriers. This may be due 
to the fact that the local environment (average number 
of kinks on the step edge) can change significantly from 
point to point along the step and also from measurement 
to measurement (as is shown quite eloquently in Fig.Qfor 
the two sets of measurements) as the equilibration time 
on larger length-scales can be very large in comparison 
to the measurement times. (See a more thorough dis- 
cussion of these experimental difficulties in M. Giesen's 
review Ref. H^|.) These measurements and especially the 
prefactors may be sensitive to these details of the system 
configuration, while other parameters seem to be more 
robust such as the dynamical exponents. 

In summary, using VT-STM we measured the funda- 
mental step edge effective energy barrier characterizing 
the transport properties of the Pb(lll) surface steps. We 
mimicked this system using KMC simulations with re- 
alistic energy barriers in the framework of the nearest 
neighbor lattice gas model. The experimental and the 
simulational data agree reasonably well for low tempera- 
ture so that we can deduce an energy barrier that shows 
the right trend for the barriers measured previously by 
other groups for Pt(lll),2i Ag(lll)^2i and Cu(lll)^ 

The use of MC simulations apparently helps greatly to 
sort out the right energy barriers and attempt frequen- 
cies from the strongly scattering experimental findings 
and places these findings on firm footings as it makes the 
trends in the data more clear. However, the comparison 
of the measurement and the MC also reveals the range 
of applicability of simple models, like the one used in 
this study with nearest neighbor interactions and (short) 
nearest neighbor jumps. While it applies to low temper- 
atures, the higher temperature parameter range appar- 
ently requires more complex models with possibly longer- 
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FIG. 8: (Color online) Three microstates and the two pro- 
cesses A — » B — > C and A — » C discussed in the text to show 
how detailed balance is satisfied for the bond-breaking energy 
barriers and fails when simply using the NN EAM barriers. 



range interactions and/or longer adatom jumps or other 
complicated atomic moves. 



APPENDIX: Detailed-balance Considerations 

To describe equilibrium thermodynamic properties ac- 
curately, Monte Carlo algorithms should satisfy detailed 
balance, which insures that the equilibrium distribution 
is a steady-state solution. Likewise, detailed balance 
is deemed crucial in simulations of systems near equi- 
librium, such as weakly driven systems. (For systems 
far from equilibrium, such as kinetically-limitcd growth, 
algorithms that violate detailed balance by excluding 
highly improbable events can still give satisfactory re- 
sults.) In this appendix we illustrate explicitly that ei- 
ther of the break schemes satisfies detailed balance. 

We consider two different ways of performing a sim- 
ple hopping process as depicted in Fig. Q3J In process 
A — > B — > C an adatom hops from site A to B and then 
C while in A — > C it hops directly from A to C. In both 
cases all the other adatoms remain fixed. A, B, and C 
can also serve to denote the microstates with the corre- 
sponding adatom configurations and energies as in Fig. 
El The energy difference of two states A.E& is the differ- 
ence between the energy barriers of the forward (A — > B) 
and backward (B — > A) moves between the states (e.g. 
AE^ B = E% 23 - El 24 in the EAM barrier notation or 

AE^ B = E ( ° ] - E^ in the brcak-3 notation) for the 
A — > B move. With a proper scheme for approximating 
the barriers, detailed balance is automatically satisfied. 
Following through this example for A — > B and B — > C 
compared to A — ► C, the two paths lead to the same en- 
ergy difference, Etd + ^Enn, between microstates A and 
C. However, doing the same with the EAM barriers the 
first and the second processes give different energy differ- 
ence between A and C: from Table |U we see that the sum 
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FIG. 9: (Color online) Energy landscape for A <-» B <-> C 
and A <-> C processes in Fig. The processes at the bottom 
and the top are indicated by the occupied NN sites around the 
hop (while in the parentheses we indicate the occupied sites 
relevant to the brake-3 model). The numbers next to the 
arrows are the EAM energy barriers in meV (in parentheses 
the corresponding break-3 barriers) taken from column four 
in Table Q (and column two for the break-3 barriers). 

-32 + 225 = 193 meV (for A -> B and B -> C processes) 
differs significantly from 235 meV (for A — > C), the dif- 
ference being 42 meV or about 20%. Thus, the NN EAM 
barriers do not satisfy detailed balance, as mentioned ear- 
lier. (When NNN's are appropriately taken into account, 
we have confirmed that detailed balance is precisely satis- 
fied. Thus, the energy landscape in EAM is well-defined. 
The subtleties arise when reducing the parameter space 
to a computationally tractable size.) 

On very general grounds, bond breaking schemes sat- 
isfy detailed balance: the energy difference between 
states i and j: 

EV j ~ = (»i - «j)£nn (9) 

where n% and rij are the number of occupied nearest- 
neighbor sites of sites i and j, respectively. Eq. © 
is a straightforward consequence of the break-5 scheme. 
Moreover, it is also valid for the break-3 case — in which 
the bridge sites 2 and 2' are neglected in m and rij — since 
these two sites are neighbors of both sites i and j; hence, 
they do not contribute to the difference of ni and rij . 

The energy difference in the situation when sites i and 
j are not neighbors but farther apart is still correctly 
given by Eq. (jj|J): the path between the initial and fi- 
nal positions can be decomposed into a (non-unique) se- 
quence of NN hops. Then the energy of each interme- 
diate position is both added and subtracted, cancelling 
out, so that difference of the n's of the initial and the 
final positions survive. Finally, since any configuration 
can be obtained from any other by moving adatoms one 
by one, every configuration has a uniquely defined energy. 
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Hence, the energy barriers form a "consistent" set. Writ- 
ing the rate of move i — ► j as = vjj exp(— fiE 1 ^ 3 ), 

where E % ^ 3 is defined in Eq. (jHJ and vd is the Debyc 
frequency, we see immediately that detailed balance is 
satisfied because the rates in the ratio 



') 



-0(n,i— ti 3 )Bnn 



(10) 



J 1 



are based on these energy levels of the configurations. 
This consistency due to the uniquely defined energies of 
the states characteristic of the bond-breaking schemes 



and also of the Kawasaki scheme is missing in the NN 
set of EAM numbers. 
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